Thunderstrom Analysis

Imports

library(evd); 
library(evdbayes);
library(coda);
library(ismev);

Lade nötiges Paket: mgcv
Lade nötiges Paket: nlme
This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.

library(xts);

Lade nötiges Paket: zoo

Attache Paket: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading the Data

load("../data/CAPE_Minder_Rychener_Malsot.RData")
load("../data/NINO34.RData")
load("../data/SRH_Minder_Rychener_Malsot.RData")

Generate PROD

prod = sqrt(cape)*srh

** Create Time Series Objects **

dates <- seq.Date(as.Date("1979-1-1"), as.Date("2015-12-31"), by=1)
feb29ix <- format(as.Date(dates), "%m-%d") == "02-29"
dates <- dates[!feb29ix]
prod_ts = xts(prod, order.by = rep(dates, each=8))

Beginning the Analysis

month_names = c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
get_monthly = function(x) {
  output = list()
  len = nrow(x)
  
  # Get month of element
  month = time(x)
  month = gsub("....-", "", month)
  month = gsub("-..", "", month)
  monthlist = unique(month)
  for (i in 1:12) {
    output[[i]] = x[month == monthlist[i],]
  }
  names(output) = month_names
  return(output)
}
monthly_max = get_monthly(apply.monthly(prod_ts, max))
r = 2
r_monthly_max = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:r]))

# get the monthly maxima of prod
m1 = as.data.frame(apply.monthly(prod_ts, max))$V1;
# produce the gumbel qq plot
gumbel_qq = function(x, title="") qqplot(qgumbel(c(1:length(x))/(length(x)+1)),
                                         x,
                                         main = paste("Gumbell Q-Q Plot", title),
                                         xlab = "Theoretical Quantiles", 
                                         ylab = "Sample Quantiles") ;
gumbel_qq(m1)

#qqplot(qgumbel(c(1:length(m1))/(length(m1)+1)),m1,main = "Gumbell Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles") ;

The qq plot doesn’t fit very well, especially in the lower tail. This is likely due to seasonal dependence.

Fitting GEV to the entire data

# fit gevd with MLE and produce diagnostic plots
fitmax.MLE<-fgev(m1,method="Nelder-Mead")
par(mfrow=c(2,2))
fitmax.MLE

Call: fgev(x = m1, method = "Nelder-Mead") 
Deviance: 9272.599 

Estimates
      loc      scale      shape  
7.519e+03  7.013e+03  4.757e-03  

Standard Errors
      loc      scale      shape  
421.90201  331.43749    0.05347  

Optimization Information
  Convergence: successful 
  Function Evaluations: 171 

plot(fitmax.MLE)

Poor fit, probably because the distribution isn’t stationary. This is especially visible in the Probability plot, in which the confidence band is surpassed, indicating a poor fit.

# fit gevd with Bayesian Techniques
# use the previous outputs (rounded) as initial values (use a different shape)
init<-c(1.6e4,4e3,0.1)
# arbitrary priors
mat <- diag(c(10000,10000,100)) 
pn <- prior.norm(mean=c(0,0,0),cov=mat)
# proposal standard deviation (found by trial and error to get 20%<acceptance rate<40%)
psd<-c(500,0.03,0.02)
# draw 3k samples from markov chain
nit = 30000
thinning = 300
post<-posterior(nit, init=init,prior=pn,lh="gev",data=m1,psd=psd)
# diagnostic plots
MCMC<-mcmc(post[1:nit %% thinning == 0, ], thin=300) 
plot(MCMC) 

attr(mcmc(post),"ar")

            mu sigma   xi total
acc.rates 0.24  0.71 0.70  0.55
ext.rates 0.00  0.00 0.01  0.00

#MCMC_stab <- mcmc(post, thin=50, start=1000)
acf(mcmc(post[1:nit %% thinning == 0, ]))

We observe that there seem to be no substantial issues from the autocorrelation plots.

apply(mcmc(post[1:nit %% thinning == 0, ]),2,mean)

           mu         sigma            xi 
  236.0524535 11181.3806761    -0.1897774 

apply(mcmc(post[1:nit %% thinning == 0, ]),2,sd)

          mu        sigma           xi 
 98.60894002 624.46021143   0.03409059 

** Fit with MLE for months separately**

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
monthly_fits = list()
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  
  if (i %in% error_cases) {
    monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  }
  else {
    monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead")
  }
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9

optimization may not have succeeded

[1] 10
[1] 11
[1] 12

names(monthly_fits) = names(monthly_max)

gumbel_qq(data.frame(monthly_max[[9]])[,1], "September")

gumbel_qq(data.frame(monthly_max[[12]])[,1], "December")

Fit R largest order statistics

largest_10 = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:4]))
largest_2_fit = lapply(largest_10, function(x) rlarg.fit(x[ , 1:2]))

$conv
[1] 0

$nllh
[1] 381.7301

$mle
[1] 172.4434267  63.1813324  -0.8246329

$se
[1] 10.3518367  6.2427031  0.1133745

$conv
[1] 0

$nllh
[1] 366.2861

$mle
[1] 155.5451810  51.4464515  -0.7267874

$se
[1] 8.3567636 4.6071290 0.1019313

$conv
[1] 0

$nllh
[1] 348.4217

$mle
[1] 195.3919929  40.8441452  -0.8324715

$se
[1] 6.33288212 4.00545299 0.08820541

$conv
[1] 0

$nllh
[1] 375.033

$mle
[1] 162.9292046  57.8831308  -0.7755148

$se
[1] 9.282767 5.432561 0.099956

$conv
[1] 0

$nllh
[1] 380.4256

$mle
[1] 157.2651459  61.7442705  -0.6554637

$se
[1] 11.5888901  5.2748342  0.1508781

$conv
[1] 0

$nllh
[1] 364.4628

$mle
[1] 104.8006071  51.1390400  -0.1253138

$se
[1] 9.153534 4.541184 0.147210

$conv
[1] 0

$nllh
[1] 386.7889

$mle
[1] 128.6412577  66.3424033  -0.4883028

$se
[1] 11.877915  5.260236  0.139483

$conv
[1] 0

$nllh
[1] 380.7543

$mle
[1] 130.1181620  62.4551520  -0.3960438

$se
[1] 10.9040496  4.7896609  0.1237002

$conv
[1] 0

$nllh
[1] 361.8412

$mle
[1] 111.0849233  49.8496920  -0.2990619

$se
[1] 8.2812262 3.8324125 0.1013576

$conv
[1] 0

$nllh
[1] 377.358

$mle
[1] 75.1622227 56.9409902  0.2522703

$se
[1] 10.0789585  7.2986557  0.1984848

$conv
[1] 0

$nllh
[1] 376.1717

$mle
[1] 106.7510273  60.2707332  -0.1192325

$se
[1] 10.4205289  5.3759827  0.1333882

$conv
[1] 0

$nllh
[1] 377.2483

$mle
[1] 126.4665127  59.0621082  -0.2939312

$se
[1] 9.9836743 4.7019303 0.1170011

lapply(largest_2_fit, rlarg.diag)

Fehler in get(".xts_chob", .plotxtsEnv) : 
  Objekt '.xts_chob' nicht gefunden

get_se = function(x, ix) {
  if (is.null(x$std.err)) 0
  else x$std.err[ix]
}
mle_loc = unlist(lapply(monthly_fits, function(x) x$estimate[1]))
mle_loc_se = unlist(lapply(monthly_fits, get_se, 1))
mle_scale = unlist(lapply(monthly_fits, function(x) x$estimate[2]))
mle_scale_se = unlist(lapply(monthly_fits, get_se, 2))
mle_shape = unlist(lapply(monthly_fits, function(x) x$estimate[3]))
mle_shape_se = unlist(lapply(monthly_fits, get_se, 3))

plot_w_err = function(x, y, se, title = NULL) {
  max_ix = which.max(y)
  min_ix = which.min(y)
  plot(x, y,
       ylim = c(y[min_ix] - se[min_ix], y[max_ix] + se[max_ix]),
       main = title)
  arrows(x,y-se,x,y+se, code=3, length=0.02, angle = 90)
}
plot_w_err(1:12, mle_loc, mle_loc_se, "Location Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_scale, mle_scale_se, "Scale Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_shape, mle_shape_se, "Shape Parameter as Estimated with Likelihood")

** Fit with Bayesian Methods for months separately**

print(acceptance_rates)

$jan
   mu sigma    xi total 
 0.31  0.32  0.33  0.32 

$feb
   mu sigma    xi total 
 0.24  0.30  0.33  0.29 

$mar
   mu sigma    xi total 
 0.24  0.32  0.20  0.25 

$apr
   mu sigma    xi total 
 0.23  0.35  0.28  0.29 

$may
   mu sigma    xi total 
 0.24  0.29  0.21  0.25 

$jun
   mu sigma    xi total 
 0.24  0.32  0.22  0.26 

$jul
   mu sigma    xi total 
 0.23  0.20  0.29  0.24 

$aug
   mu sigma    xi total 
 0.24  0.21  0.25  0.23 

$sep
   mu sigma    xi total 
 0.23  0.27  0.20  0.23 

$oct
   mu sigma    xi total 
 0.24  0.38  0.23  0.28 

$nov
   mu sigma    xi total 
 0.24  0.34  0.34  0.31 

$dec
   mu sigma    xi total 
 0.34  0.33  0.30  0.33 

TODO-> R largest fit

PART 2 First, we check if the location parameter depends on time using a likelihood ratio test

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  
  fit_const = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             nsloc = trend,
                             std.err = FALSE)
  
  ratios[[i]] = fit_const$dev-fit_dependant$dev 
}

[1] 1

optimization may not have succeeded

[1] 2

optimization may not have succeeded

[1] 3
[1] 4

optimization may not have succeeded

[1] 5

optimization may not have succeeded

[1] 6
[1] 7
[1] 8
[1] 9

optimization may not have succeeded

[1] 10
[1] 11

optimization may not have succeeded

[1] 12

names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for time independance, Bonferroni multiple Testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Now, let’s check for independance from ENSO

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
# split nino data into months
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  trend = n[i,]
  trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
  fit_const = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             nsloc = trend,
                             std.err = FALSE)
  
  ratios[[i]] = fit_const$dev-fit_dependant$dev 
}

[1] 1

optimization may not have succeeded

[1] 2

optimization may not have succeeded

[1] 3
[1] 4

optimization may not have succeeded

[1] 5
[1] 6

optimization may not have succeeded

[1] 7
[1] 8

optimization may not have succeeded

[1] 9

optimization may not have succeeded

[1] 10
[1] 11
[1] 12

names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for independance from ENSO, Bonferroni multiple testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Another method is the chi plot:

# plot the chi plot for dependance analysis
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
for (i in 1:length(monthly_max)) {
  nino = n[i,]
  m_data=as.data.frame(monthly_max[[i]])$V1
  dat.m1_month = cbind(m_data,trend);
  dat.m1_nino = cbind(m_data,nino);
  par(mfrow=c(2,2))
  chiplot(dat.m1_month,main1 = "Chi Plot Time",main2 = "Chi Bar Plot Time");
  chiplot(dat.m1_nino,main1 = "Chi Plot ENSO",main2 = "Chi Bar Plot ENSO");
}

PART 3 We will now analyse temporal clustering of extremes. For this, we will use the exiplot function from the evd library.

# define a function for getting the extremal indices for each month for a given threshold
monthly_eindex <- function(data, threshold_p, r=0){
  ei = list()
  for (i in 1:length(data)) {
    threshold = quantile(as.data.frame(data[[i]])$V1, threshold_p)
    ei[[i]]=exi(as.data.frame(data[[i]])$V1, threshold, r)
  }
  names(ei) = names(data)
  return(ei)
}
ei = monthly_eindex(get_monthly(prod_ts), 0.95)
plot(unlist(ei), main="Extremal Index by Month, 95%-Quantile as Threshold", xlab="Month", ylab="Extremal index")

We observe that the extremal index is _0.25-0.45, we can therefore conclude that we have strong dependance of extremes, with the limiting mean cluster size being roughly from 2 to 4. The clustering has no effect for estimaters based on the (monthly) maximum, but the r largest estimator is influenced by it.

PART 4 First, let’s estimate the 10 year return level using point process approach

monthly_fits_pp = list()
monthly_data = get_monthly(prod_ts)
error_cases = c(5,6,7,8,9,10,11,12)
month_days = c(31,28,31,30,31,30,31,31,30,31,30,31)
for (i in 1:length(monthly_max)) {
  print(i)
  threshold = quantile(as.data.frame(monthly_data[[i]])$V1, 0.95)
  
  if (i %in% error_cases) {
    monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
                             threshold = threshold,
                             model="pp",
                             npp = month_days[i]*8,
                             cmax = TRUE,
                             r = 1,
                             std.err = FALSE,
                             method = "Nelder-Mead")
  }
  else {
    monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
                             threshold = threshold,
                             model="pp",
                             npp = month_days[i]*8,
                             cmax = TRUE,
                             r = 1,
                             method = "Nelder-Mead")
  }
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12

names(monthly_fits_pp) = names(monthly_data)
for(i in 1:12){
  par(mfrow=c(2,2)) 
  plot(monthly_fits_pp[[i]])
}

The fit in february has completely failed, and the others are not very good either

We will still estimate the return levels:

return_level = function(x,period=20){
  p = 1/period
  loc = x$estimate[[1]]
  scale = x$estimate[[2]]
  shape = x$estimate[[3]]
  level = loc + scale*(((-log(1-p))^-shape-1)/shape)
  return(level)
}
return_level_20 = lapply(monthly_fits, return_level) # 20 for testing
return_level_100 = lapply(monthly_fits, return_level, period=100)
return_level_50 = lapply(monthly_fits, return_level, period=50)
plot(unlist(return_level_100),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

plot(unlist(return_level_50),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

TODO: estimat with mcmc Assuming that we have the posterior densities of the markov chains, call theis function to plot return level histograms

return_level_mcmc = function(posterior,period=20){
  u = mc.quant(posterior,p=1-1/period,lh="gev")
  label_mcmc_rl = sprintf("%s Year Return Level",period)
  hist(u,nclass=20,prob=T,xlab=label_mcmc_rl, main = "Return Level Histogram")
}
lapply(monthly_bayes_fit, function(x) return_level_mcmc(x$posterior))

$jan
$breaks
 [1]     0  5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000 75000 80000 85000 90000

$counts
 [1] 14646 12954  1775   383   131    58    15    17     7     4     3     3     0     1
[15]     1     2     0     1

$density
 [1] 9.763675e-05 8.635712e-05 1.183294e-05 2.553248e-06 8.733042e-07 3.866538e-07
 [7] 9.999667e-08 1.133296e-07 4.666511e-08 2.666578e-08 1.999933e-08 1.999933e-08
[13] 0.000000e+00 6.666444e-09 6.666444e-09 1.333289e-08 0.000000e+00 6.666444e-09

$mids
 [1]  2500  7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
[15] 72500 77500 82500 87500

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$feb
$breaks
 [1]     0  5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000

$counts
 [1]     1 18276 10221  1129   242    74    32    14     5     2     2     1     2

$density
 [1] 6.666444e-09 1.218359e-04 6.813773e-05 7.526416e-06 1.613280e-06 4.933169e-07
 [7] 2.133262e-07 9.333022e-08 3.333222e-08 1.333289e-08 1.333289e-08 6.666444e-09
[13] 1.333289e-08

$mids
 [1]  2500  7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$mar
$breaks
 [1]  4000  6000  8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000 46000 48000 50000 52000 54000

$counts
 [1]     1     1     0     0     0     0   103  1680  7830 11503  5998  1980   593   178
[15]    66    35    16    11     2     2     1     0     0     0     1

$density
 [1] 1.666611e-08 1.666611e-08 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
 [7] 1.716609e-06 2.799907e-05 1.304957e-04 1.917103e-04 9.996333e-05 3.299890e-05
[13] 9.883004e-06 2.966568e-06 1.099963e-06 5.833139e-07 2.666578e-07 1.833272e-07
[19] 3.333222e-08 3.333222e-08 1.666611e-08 0.000000e+00 0.000000e+00 0.000000e+00
[25] 1.666611e-08

$mids
 [1]  5000  7000  9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000 45000 47000 49000 51000 53000

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$apr
$breaks
 [1]     0  5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000

$counts
 [1]     1     3     3     3    88 20707  8892   270    21     7     1     3     1     1

$density
 [1] 6.666444e-09 1.999933e-08 1.999933e-08 1.999933e-08 5.866471e-07 1.380421e-04
 [7] 5.927802e-05 1.799940e-06 1.399953e-07 4.666511e-08 6.666444e-09 1.999933e-08
[13] 6.666444e-09 6.666444e-09

$mids
 [1]  2500  7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$may
$breaks
 [1]     0  5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000 75000 80000 85000

$counts
 [1]     1     3     3     3     1     1  1979 24216  3640   125    15     7     2     2
[15]     2     0     1

$density
 [1] 6.666444e-09 1.999933e-08 1.999933e-08 1.999933e-08 6.666444e-09 6.666444e-09
 [7] 1.319289e-05 1.614346e-04 2.426586e-05 8.333056e-07 9.999667e-08 4.666511e-08
[13] 1.333289e-08 1.333289e-08 1.333289e-08 0.000000e+00 6.666444e-09

$mids
 [1]  2500  7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
[15] 72500 77500 82500

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$jun
$breaks
 [1]     0  5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000 75000 80000 85000

$counts
 [1]     1     3     3     3     1     1  2192 22381  5201   185    17     6     2     2
[15]     2     0     1

$density
 [1] 6.666444e-09 1.999933e-08 1.999933e-08 1.999933e-08 6.666444e-09 6.666444e-09
 [7] 1.461285e-05 1.492017e-04 3.467218e-05 1.233292e-06 1.133296e-07 3.999867e-08
[13] 1.333289e-08 1.333289e-08 1.333289e-08 0.000000e+00 6.666444e-09

$mids
 [1]  2500  7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
[15] 72500 77500 82500

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$jul
$breaks
 [1]  4000  6000  8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000

$counts
 [1]     2     3     1     3     2     0     1     0   522 13995 14730   663    54    15
[15]     6     2     1     0     1

$density
 [1] 3.333222e-08 4.999833e-08 1.666611e-08 4.999833e-08 3.333222e-08 0.000000e+00
 [7] 1.666611e-08 0.000000e+00 8.699710e-06 2.332422e-04 2.454918e-04 1.104963e-05
[13] 8.999700e-07 2.499917e-07 9.999667e-08 3.333222e-08 1.666611e-08 0.000000e+00
[19] 1.666611e-08

$mids
 [1]  5000  7000  9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$aug
$breaks
 [1]  4000  6000  8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000

$counts
 [1]     1     1     0     0     2     0     0     6   216  3344 12190 10653  2914   518
[15]   115    20    14     2     2     3

$density
 [1] 1.666611e-08 1.666611e-08 0.000000e+00 0.000000e+00 3.333222e-08 0.000000e+00
 [7] 0.000000e+00 9.999667e-08 3.599880e-06 5.573148e-05 2.031599e-04 1.775441e-04
[13] 4.856505e-05 8.633046e-06 1.916603e-06 3.333222e-07 2.333256e-07 3.333222e-08
[19] 3.333222e-08 4.999833e-08

$mids
 [1]  5000  7000  9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$sep
$breaks
 [1]     0  5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000

$counts
 [1]     1     3     3    15 23255  6680    32     4     0     2     2     3     1

$density
 [1] 6.666444e-09 1.999933e-08 1.999933e-08 9.999667e-08 1.550282e-04 4.453185e-05
 [7] 2.133262e-07 2.666578e-08 0.000000e+00 1.333289e-08 1.333289e-08 1.999933e-08
[13] 6.666444e-09

$mids
 [1]  2500  7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$oct
$breaks
 [1]  4000  6000  8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000 46000 48000

$counts
 [1]    1    1    0    0    0   30  938 5430 9239 8150 3932 1437  485  204   70   44   20
[18]    9    3    3    1    4

$density
 [1] 1.666611e-08 1.666611e-08 0.000000e+00 0.000000e+00 0.000000e+00 4.999833e-07
 [7] 1.563281e-05 9.049698e-05 1.539782e-04 1.358288e-04 6.553115e-05 2.394920e-05
[13] 8.083064e-06 3.399887e-06 1.166628e-06 7.333089e-07 3.333222e-07 1.499950e-07
[19] 4.999833e-08 4.999833e-08 1.666611e-08 6.666444e-08

$mids
 [1]  5000  7000  9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000 45000 47000

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$nov
$breaks
 [1]  4000  6000  8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000

$counts
 [1]     1   102  5217 14247  7181  2090   689   269   101    50    19    14     8     1
[15]     2     5     2     2     0     1

$density
 [1] 1.666611e-08 1.699943e-06 8.694710e-05 2.374421e-04 1.196793e-04 3.483217e-05
 [7] 1.148295e-05 4.483184e-06 1.683277e-06 8.333056e-07 3.166561e-07 2.333256e-07
[13] 1.333289e-07 1.666611e-08 3.333222e-08 8.333056e-08 3.333222e-08 3.333222e-08
[19] 0.000000e+00 1.666611e-08

$mids
 [1]  5000  7000  9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

$dec
$breaks
 [1]     0  5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000

$counts
 [1] 23735  5642   470   103    16    18     9     1     2     0     0     0     3     2

$density
 [1] 1.582281e-04 3.761208e-05 3.133229e-06 6.866438e-07 1.066631e-07 1.199960e-07
 [7] 5.999800e-08 6.666444e-09 1.333289e-08 0.000000e+00 0.000000e+00 0.000000e+00
[13] 1.999933e-08 1.333289e-08

$mids
 [1]  2500  7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500

$xname
[1] "u"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"